
#############################################################################
# Public Support for Democracy and Three Patterns of Democratic Backsliding
# Author: Jeongho Choi and Byung-deuk Woo
# Last Updated: Feb 12nd, 2022
# Purpose: Merge Data & Run the Model
# Version: R version 4.1.1 -- "Kick Things"
############################################################################

# 0. Set ================================================
## 0.1 Setting the Working Directory=====================
rm(list=ls())
directory<-getwd()
setwd(directory)
memory.limit(size = 56000)

## 0.2 Packages ========================================
library(rio)
library(dplyr)
library(tidyverse)
library(stargazer)
library(BayesFactor)
library(MCMCpack)
library(MCMCglmm)
library(BMS)
library(BAS)
library(MASS)
library(arm)
library(tidyr)
library(loo)
library(ggplot2)
library(arm)
library(bayesplot)
library(brms)
library(fs)
library(dotwhisker)
library(ggplot2)

## 0.3  Importing Data ======================================================
#unzip(zipfile = "Final_Democratic Backsliding.zip", exdir = "./data_zip")
#data_dir<-"data_zip"

csv_files<-fs::dir_ls(directory, regexp = "\\.csv$")
data<-csv_files %>% map_dfr(read_csv)
unlink(data_dir, recursive = TRUE)

#Excluded closed democracy
data<-data[!(is.na(data$regime_type)==TRUE),]
data<-data[!(data$regime_type==0),]

#Split_regime
data_dem<-data[(data$regime_type==2|data$regime_type==3),]
data_auto<-data[(data$regime_type==1),]


## 1. Bayesian Regression ###################################################
## 1.1 Claassen Model =======================================================
exec_all_claassen_nofix =lm(vdem_horacc ~ vdem_horacc_lag1 + vdem_horacc_lag2
                      +supdem_lag1
                      +log(gdp_pc_lag1) # GDP per Capita
                      +gdp_pc_growth_lag1 # GDP annual growth
                      +natural_resource_lag1 # Natural Resource
                      +dem_region_lag1, # Mean value of Liberal democracy by Region
                      data)

elec_all_claassen_nofix =lm(vdem_frefair ~ vdem_frefair_lag1 + vdem_frefair_lag2
                      +supdem_lag1
                      +log(gdp_pc_lag1) # GDP per Capita
                      +gdp_pc_growth_lag1 # GDP annual growth
                      +natural_resource_lag1 # Natural Resource
                      +dem_region_lag1, # Mean value of Liberal democracy by Region
                      data)

bayes_exec_all_claassen_nofix<-MCMCregress(formula =exec_all_claassen_nofix, b0=c(0, #Intercept
                                                                      0, #v_dem lag1
                                                                      0, #v_dem lag2
                                                                      0, #Supdem
                                                                      0, # GDPpc
                                                                      0, # GDPpc growth
                                                                      0, #Natural Resource
                                                                      0),#Dem region
                                     B0=0.01, data=data, marginal.likelihood = "Laplace",
                                     mcmc = 10000, burnin = 2000, seed = 12345)

bayes_elec_all_claassen_nofix<-MCMCregress(formula =elec_all_claassen_nofix, b0=c(0, #Intercept
                                                                      0, #v_dem lag1
                                                                      0, #v_dem lag2
                                                                      0, #Supdem
                                                                      0, # GDPpc
                                                                      0, # GDPpc growth
                                                                      0, #Natural Resource
                                                                      0),#Dem region
                                     B0=0.01, data=data, marginal.likelihood = "Laplace",
                                     mcmc = 10000, burnin = 2000, seed = 12345)

#In democracy
exec_dem_claassen_nofix =lm(vdem_horacc ~ vdem_horacc_lag1 + vdem_horacc_lag2
                      +supdem_lag1
                      +log(gdp_pc_lag1) # GDP per Capita
                      +gdp_pc_growth_lag1 # GDP annual growth
                      +natural_resource_lag1 # Natural Resource
                      +dem_region_lag1, # Mean value of Liberal democracy by Region
                      data_dem)
elec_dem_claassen_nofix =lm(vdem_frefair ~ vdem_frefair_lag1 + vdem_frefair_lag2
                            +supdem_lag1
                            +log(gdp_pc_lag1) # GDP per Capita
                            +gdp_pc_growth_lag1 # GDP annual growth
                            +natural_resource_lag1 # Natural Resource
                            +dem_region_lag1, # Mean value of Liberal democracy by Region
                      data_dem)

bayes_exec_dem_claassen_nofix<-MCMCregress(formula =exec_dem_claassen_nofix, b0=c(0, #Intercept
                                                                      0, #v_dem lag1
                                                                      0, #v_dem lag2
                                                                      0, #Supdem
                                                                      0, # GDPpc
                                                                      0, # GDPpc growth
                                                                      0, #Natural Resource
                                                                      0),#Dem region
                                     B0=0.01, data=data_dem, marginal.likelihood = "Laplace",
                                     mcmc = 10000, burnin = 2000, seed = 12345)

bayes_elec_dem_claassen_nofix<-MCMCregress(formula =elec_dem_claassen_nofix, b0=c(0, #Intercept
                                                                      0, #v_dem lag1
                                                                      0, #v_dem lag2
                                                                      0, #Supdem
                                                                      0, # GDPpc
                                                                      0, # GDPpc growth
                                                                      0, #Natural Resource
                                                                      0),#Dem region
                                     B0=0.01, data=data_dem, marginal.likelihood = "Laplace",
                                     mcmc = 10000, burnin = 2000, seed = 12345)

#In autocracy
exec_auto_claassen_nofix =lm(vdem_horacc ~ vdem_horacc_lag1 + vdem_horacc_lag2
                       +supdem_lag1
                       +log(gdp_pc_lag1) # GDP per Capita
                       +gdp_pc_growth_lag1 # GDP annual growth
                       +natural_resource_lag1 # Natural Resource
                       +dem_region_lag1, # Mean value of Liberal democracy by Region
                       data_auto)
elec_auto_claassen_nofix =lm(vdem_frefair ~ vdem_frefair_lag1 + vdem_frefair_lag2
                             +supdem_lag1
                             +log(gdp_pc_lag1) # GDP per Capita
                             +gdp_pc_growth_lag1 # GDP annual growth
                             +natural_resource_lag1 # Natural Resource
                             +dem_region_lag1, # Mean value of Liberal democracy by Region
                       data_auto)

bayes_exec_auto_claassen_nofix<-MCMCregress(formula =exec_auto_claassen_nofix, b0=c(0, #Intercept
                                                                        0, #v_dem lag1
                                                                        0, #v_dem lag2
                                                                        0, #Supdem
                                                                        0, # GDPpc
                                                                        0, # GDPpc growth
                                                                        0, #Natural Resource
                                                                        0),#Dem region
                                      B0=0.01, data=data_auto, marginal.likelihood = "Laplace",
                                      mcmc = 10000, burnin = 2000, seed = 12345)

bayes_elec_auto_claassen_nofix<-MCMCregress(formula =elec_auto_claassen_nofix, b0=c(0, #Intercept
                                                                        0, #v_dem lag1
                                                                        0, #v_dem lag2
                                                                        0, #Supdem
                                                                        0, # GDPpc
                                                                        0, # GDPpc growth
                                                                        0, #Natural Resource
                                                                        0),#Dem region
                                      B0=0.01, data=data_auto, marginal.likelihood = "Laplace",
                                      mcmc = 10000, burnin = 2000, seed = 12345)

## 2.Result ###################################################################
## 2.1 Result Table ==========================================================
#Table 1
summary(bayes_exec_all_claassen_nofix ) #Negative
summary(bayes_elec_all_claassen_nofix ) #Negative
a<-summary(bayes_exec_all_claassen_nofix)$quantiles[,c(1,3,5)]
b<-summary(bayes_elec_all_claassen_nofix)$quantiles[,c(1,3,5)]
stargazer(a,b, type = "text",digits = 2)

summary(bayes_exec_dem_claassen_nofix ) #Negative
summary(bayes_elec_dem_claassen_nofix ) #Negative
a<-summary(bayes_exec_dem_claassen_nofix)$quantiles[,c(1,3,5)]
b<-summary(bayes_elec_dem_claassen_nofix)$quantiles[,c(1,3,5)]
stargazer(a,b, type = "text",digits = 2)

summary(bayes_exec_auto_claassen_nofix) #Negative
summary(bayes_elec_auto_claassen_nofix) #Null
a<-summary(bayes_exec_auto_claassen_nofix)$quantiles[,c(1,3,5)]
b<-summary(bayes_elec_auto_claassen_nofix)$quantiles[,c(1,3,5)]
stargazer(a,b, type = "text",digits = 2)

## 2.2 Visualization ============================================================================
#Figure 1
plot(density(bayes_exec_all_claassen_nofix[ , 4]),xlim = c(-0.003, 0), ylim=c(0,7000),
     lwd = 2, lty =1,  col = "black", xlab = "Coefficient of public support for democracy",
     main = "")
lines(density( bayes_exec_dem_claassen_nofix [ , 4]), lwd = 2, lty=1, col = "blue")
lines(density(bayes_exec_auto_claassen_nofix [ , 4]), lwd = 2, lty=1, col = "red")
abline(v=0,lwd=3, lty=2)
legend("topright", lty=1, lwd=2,
       legend=c("General", "In democracy", "In non-democracy"), col=c("black", "blue", "red"))

#Figure 2
plot(density(bayes_elec_all_claassen_nofix[ , 4]),xlim = c(-0.001, 0.0001), ylim=c(0,15000),
     lwd = 2, lty =1,  col = "black", xlab = "Coefficient of public support for democracy",
     main = "")
lines(density( bayes_elec_dem_claassen_nofix [ , 4]), lwd = 2, lty=1, col = "blue")
lines(density( bayes_elec_auto_claassen_nofix [ , 4]), lwd = 2, lty=1, col = "red")
abline(v=0,lwd=3, lty=2)
legend("topright", lty=1, lwd=2,
       legend=c("General", "In democracy", "In non-democracy"), col=c("black", "blue", "red"))

## 3 Robustness test ###########################################################
## 3.1 Convergence diagnostics =================================================
#Traceplots
par(mfrow=c(3,3))
traceplot(bayes_exec_all_claassen_nofix)
traceplot(bayes_elec_all_claassen_nofix)

#Running mean
cumuplot(bayes_exec_all_claassen_nofix, probs=c(.5))
cumuplot(bayes_elec_all_claassen_nofix, probs=c(.5))

#Geweke plot
geweke.plot(bayes_exec_all_claassen_nofix)
geweke.plot(bayes_elec_all_claassen_nofix)

#Autocorrelation
autocorr.plot(bayes_exec_all_claassen_nofix, 50)
autocorr.plot(bayes_elec_all_claassen_nofix, 50)

## 3.2 Bayes Factor Analysis ===================================================
exec_all_claassen_nofix =lm(vdem_horacc ~ vdem_horacc_lag1 + vdem_horacc_lag2
                            +supdem_lag1
                            +log(gdp_pc_lag1) # GDP per Capita
                            +gdp_pc_growth_lag1 # GDP annual growth
                            +natural_resource_lag1 # Natural Resource
                            +dem_region_lag1, # Mean value of Liberal democracy by Region
                            data)
elec_all_claassen_nofix =lm(vdem_frefair ~ vdem_frefair_lag1 + vdem_frefair_lag2
                            +supdem_lag1
                            +log(gdp_pc_lag1) # GDP per Capita
                            +gdp_pc_growth_lag1 # GDP annual growth
                            +natural_resource_lag1 # Natural Resource
                            +dem_region_lag1, # Mean value of Liberal democracy by Region
                            data)

bayes_exec_all_claassen_nofix_positive<-MCMCregress(formula =exec_all_claassen_nofix, b0=c(0, #Intercept
                                                                                  0, #v_dem lag1
                                                                                  0, #v_dem lag2
                                                                                  1, #Supdem
                                                                                  0, # GDPpc
                                                                                  0, # GDPpc growth
                                                                                  0, #Natural Resource
                                                                                  0),#Dem region
                                           B0=0.01, data=data, marginal.likelihood = "Laplace",
                                           mcmc = 10000, burnin = 2000, seed = 12345)
bayes_exec_all_claassen_nofix_negative<-MCMCregress(formula =exec_all_claassen_nofix, b0=c(0, #Intercept
                                                                                           0, #v_dem lag1
                                                                                           0, #v_dem lag2
                                                                                           -1, #Supdem
                                                                                           0, # GDPpc
                                                                                           0, # GDPpc growth
                                                                                           0, #Natural Resource
                                                                                           0),#Dem region
                                                    B0=0.01, data=data, marginal.likelihood = "Laplace",
                                                    mcmc = 10000, burnin = 2000, seed = 12345)

bayes_elec_all_claassen_nofix_positive<-MCMCregress(formula =elec_all_claassen_nofix, b0=c(0, #Intercept
                                                                                  0, #v_dem lag1
                                                                                  0, #v_dem lag2
                                                                                  1, #Supdem
                                                                                  0, # GDPpc
                                                                                  0, # GDPpc growth
                                                                                  0, #Natural Resource
                                                                                  0),#Dem region
                                           B0=0.01, data=data, marginal.likelihood = "Laplace",
                                           mcmc = 10000, burnin = 2000, seed = 12345)

bayes_elec_all_claassen_nofix_negative<-MCMCregress(formula =elec_all_claassen_nofix, b0=c(0, #Intercept
                                                                                  0, #v_dem lag1
                                                                                  0, #v_dem lag2
                                                                                  -1, #Supdem
                                                                                  0, # GDPpc
                                                                                  0, # GDPpc growth
                                                                                  0, #Natural Resource
                                                                                  0),#Dem region
                                           B0=0.01, data=data, marginal.likelihood = "Laplace",
                                           mcmc = 10000, burnin = 2000, seed = 12345)
#In democracy
exec_dem_claassen_nofix =lm(vdem_horacc ~ vdem_horacc_lag1 + vdem_horacc_lag2
                            +supdem_lag1
                            +log(gdp_pc_lag1) # GDP per Capita
                            +gdp_pc_growth_lag1 # GDP annual growth
                            +natural_resource_lag1 # Natural Resource
                            +dem_region_lag1, # Mean value of Liberal democracy by Region
                            data_dem)
elec_dem_claassen_nofix =lm(vdem_frefair ~ vdem_frefair_lag1 + vdem_frefair_lag2
                            +supdem_lag1
                            +log(gdp_pc_lag1) # GDP per Capita
                            +gdp_pc_growth_lag1 # GDP annual growth
                            +natural_resource_lag1 # Natural Resource
                            +dem_region_lag1, # Mean value of Liberal democracy by Region
                            data_dem)

bayes_exec_dem_claassen_nofix_positive<-MCMCregress(formula =exec_dem_claassen_nofix, b0=c(0, #Intercept
                                                                                  0, #v_dem lag1
                                                                                  0, #v_dem lag2
                                                                                  1, #Supdem
                                                                                  0, # GDPpc
                                                                                  0, # GDPpc growth
                                                                                  0, #Natural Resource
                                                                                  0),#Dem region
                                           B0=0.01, data=data_dem, marginal.likelihood = "Laplace",
                                           mcmc = 10000, burnin = 2000, seed = 12345)
bayes_exec_dem_claassen_nofix_negative<-MCMCregress(formula =exec_dem_claassen_nofix, b0=c(0, #Intercept
                                                                                  0, #v_dem lag1
                                                                                  0, #v_dem lag2
                                                                                  -1, #Supdem
                                                                                  0, # GDPpc
                                                                                  0, # GDPpc growth
                                                                                  0, #Natural Resource
                                                                                  0),#Dem region
                                           B0=0.01, data=data_dem, marginal.likelihood = "Laplace",
                                           mcmc = 10000, burnin = 2000, seed = 12345)

bayes_elec_dem_claassen_nofix_positive<-MCMCregress(formula =elec_dem_claassen_nofix, b0=c(0, #Intercept
                                                                                  0, #v_dem lag1
                                                                                  0, #v_dem lag2
                                                                                  1, #Supdem
                                                                                  0, # GDPpc
                                                                                  0, # GDPpc growth
                                                                                  0, #Natural Resource
                                                                                  0),#Dem region
                                           B0=0.01, data=data_dem, marginal.likelihood = "Laplace",
                                           mcmc = 10000, burnin = 2000, seed = 12345)
bayes_elec_dem_claassen_nofix_negative<-MCMCregress(formula =elec_dem_claassen_nofix, b0=c(0, #Intercept
                                                                                  0, #v_dem lag1
                                                                                  0, #v_dem lag2
                                                                                  -1, #Supdem
                                                                                  0, # GDPpc
                                                                                  0, # GDPpc growth
                                                                                  0, #Natural Resource
                                                                                  0),#Dem region
                                           B0=0.01, data=data_dem, marginal.likelihood = "Laplace",
                                           mcmc = 10000, burnin = 2000, seed = 12345)
#In autocracy
exec_auto_claassen_nofix =lm(vdem_horacc ~ vdem_horacc_lag1 + vdem_horacc_lag2
                             +supdem_lag1
                             +log(gdp_pc_lag1) # GDP per Capita
                             +gdp_pc_growth_lag1 # GDP annual growth
                             +natural_resource_lag1 # Natural Resource
                             +dem_region_lag1, # Mean value of Liberal democracy by Region
                             data_auto)
elec_auto_claassen_nofix =lm(vdem_frefair ~ vdem_frefair_lag1 + vdem_frefair_lag2
                             +supdem_lag1
                             +log(gdp_pc_lag1) # GDP per Capita
                             +gdp_pc_growth_lag1 # GDP annual growth
                             +natural_resource_lag1 # Natural Resource
                             +dem_region_lag1, # Mean value of Liberal democracy by Region
                             data_auto)

bayes_exec_auto_claassen_nofix_positive<-MCMCregress(formula =exec_auto_claassen_nofix, b0=c(0, #Intercept
                                                                                    0, #v_dem lag1
                                                                                    0, #v_dem lag2
                                                                                    1, #Supdem
                                                                                    0, # GDPpc
                                                                                    0, # GDPpc growth
                                                                                    0, #Natural Resource
                                                                                    0),#Dem region
                                            B0=0.01, data=data_auto, marginal.likelihood = "Laplace",
                                            mcmc = 10000, burnin = 2000, seed = 12345)
bayes_exec_auto_claassen_nofix_negative<-MCMCregress(formula =exec_auto_claassen_nofix, b0=c(0, #Intercept
                                                                                    0, #v_dem lag1
                                                                                    0, #v_dem lag2
                                                                                    -1, #Supdem
                                                                                    0, # GDPpc
                                                                                    0, # GDPpc growth
                                                                                    0, #Natural Resource
                                                                                    0),#Dem region
                                            B0=0.01, data=data_auto, marginal.likelihood = "Laplace",
                                            mcmc = 10000, burnin = 2000, seed = 12345)
bayes_elec_auto_claassen_nofix_positive<-MCMCregress(formula =elec_auto_claassen_nofix, b0=c(0, #Intercept
                                                                                    0, #v_dem lag1
                                                                                    0, #v_dem lag2
                                                                                    1, #Supdem
                                                                                    0, # GDPpc
                                                                                    0, # GDPpc growth
                                                                                    0, #Natural Resource
                                                                                    0),#Dem region
                                            B0=0.01, data=data_auto, marginal.likelihood = "Laplace",
                                            mcmc = 10000, burnin = 2000, seed = 12345)
bayes_elec_auto_claassen_nofix_negative<-MCMCregress(formula =elec_auto_claassen_nofix, b0=c(0, #Intercept
                                                                                    0, #v_dem lag1
                                                                                    0, #v_dem lag2
                                                                                    -1, #Supdem
                                                                                    0, # GDPpc
                                                                                    0, # GDPpc growth
                                                                                    0, #Natural Resource
                                                                                    0),#Dem region
                                            B0=0.01, data=data_auto, marginal.likelihood = "Laplace",
                                            mcmc = 10000, burnin = 2000, seed = 12345)


BayesFactor(bayes_exec_all_claassen_nofix,
            bayes_exec_all_claassen_nofix_positive,
            bayes_exec_all_claassen_nofix_negative)$BF.log.mat

BayesFactor(bayes_elec_all_claassen_nofix,
            bayes_elec_all_claassen_nofix_positive,
            bayes_elec_all_claassen_nofix_negative)$BF.log.mat

BayesFactor(bayes_exec_dem_claassen_nofix,
            bayes_exec_dem_claassen_nofix_positive,
            bayes_exec_dem_claassen_nofix_negative)$BF.log.mat

BayesFactor(bayes_elec_dem_claassen_nofix,
            bayes_elec_dem_claassen_nofix_positive,
            bayes_elec_dem_claassen_nofix_negative)$BF.log.mat

BayesFactor(bayes_exec_auto_claassen_nofix,
            bayes_exec_auto_claassen_nofix_positive,
            bayes_exec_auto_claassen_nofix_negative)$BF.log.mat

BayesFactor(bayes_elec_auto_claassen_nofix,
            bayes_elec_auto_claassen_nofix_positive,
            bayes_elec_auto_claassen_nofix_negative)$BF.log.mat

## 3.3 Prior Sensitive Test=====================================================
#Prior Sensitive Test
summary(bayes_exec_all_claassen_nofix_positive ) #Negative
summary(bayes_elec_all_claassen_nofix_positive ) #Negative
summary(bayes_exec_dem_claassen_nofix_positive ) #Negative
summary(bayes_elec_dem_claassen_nofix_positive ) #Negative
summary(bayes_exec_auto_claassen_nofix_positive) #Negative
summary(bayes_elec_auto_claassen_nofix_positive) #Null


summary(bayes_exec_all_claassen_nofix_negative ) #Negative
summary(bayes_elec_all_claassen_nofix_negative ) #Negative
summary(bayes_exec_dem_claassen_nofix_negative ) #Negative
summary(bayes_elec_dem_claassen_nofix_negative ) #Negative
summary(bayes_exec_auto_claassen_nofix_negative) #Negative
summary(bayes_elec_auto_claassen_nofix_negative) #Null

plot(density(bayes_exec_all_claassen_nofix[ , 2]), lwd = 2, lty =1,  col = "black", xlab = "beta",
     main = "Posterior Density: Constraints on the Exeuctive")
lines(density(bayes_exec_all_claassen_nofix_positive[ , 2]), lwd = 2, lty=1, col = "blue")
lines(density(bayes_exec_all_claassen_nofix_negative[ , 2]), lwd = 2, lty=1, col = "red")
legend("topright", lty=1, lwd=2, 
       legend=c("mu= 0", "mu= 10", "mu= -10"), col=c("black", "blue", "red"))

plot(density(bayes_elec_all_claassen_nofix[ , 2]), lwd = 2, lty =1,  col = "black", xlab = "beta",
     main = "Posterior Density: Clean Election")
lines(density(bayes_elec_all_claassen_nofix_positive[ , 2]), lwd = 2, lty=1, col = "blue")
lines(density(bayes_elec_all_claassen_nofix_negative[ , 2]), lwd = 2, lty=1, col = "red")
legend("topright", lty=1, lwd=2, 
       legend=c("mu= 0", "mu= 10", "mu= -10"), col=c("black", "blue", "red"))

## 3.4 Summary Statistics ======================================================
summary_table<-data[,c("vdem_horacc", "vdem_frefair",
                       "gdp_pc_lag1", "gdp_pc_growth_lag1", "natural_resource_lag1",
                       "dem_region_lag1", "regime_type")]
summary_table<-unique(summary_table)
summary_table$log_gdp_pc_lag1<-log(summary_table$gdp_pc_lag1)
stargazer(as.data.frame(summary_table),
          type = "text",digits = 2)

summary_supdem<-data$supdem

stargazer(as.data.frame(summary_supdem),
          type = "text",digits = 2)

hist(summary_supdem)